Classification of initial state granularity via 2d Fourier expansion 
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A new method to quantify fluctuations in the initial state of heavy ion collisions is presented. 
The initial state energy distribution is decomposed with a set of orthogonal basis functions which 
include both angular and radial variation. The resulting two dimensional Fourier coefficients pro- 
vide additional information about the nature of the initial state fluctuations compared to a purely 
angular decomposition. We apply this method to ensembles of initial states generated by both 
Glauber and Color Glass Condensate Monte-Carlo codes. In addition initial state configurations 
with varying amounts of fluctuations generated by a dynamic transport approach are analysed to 
test the sensitivity of the procedure. The results allow for a full characterization of the initial state 
structures that is useful to discriminate the different initial state models currently in use. 



Ultra-Relativistic nearly-ideal fluid dynamics has 
proven to be a very successful tool for modeling the bulk 
dynamics of the hot dense matter formed during a heavy 
ion collision \V^ . The major uncertainty in determining 
transport properties of the QGP, such as the shear viscos- 
ity to entropy ratio, lies in the specification of the initial 
conditions of the collision. The initial conditions have 
been mainly assumed to be smooth distributions that 
are parametrized implementations of certain physical as- 
sumptions (e.g., Glauber/CGC). Within the last 2 years 
the importance of including fluctuations in these distri- 
butions has been recognized, leading to a whole new set of 
experimental observations of higher flow coefficients and 
their correlations [THTU], On the theoretical side there 
has been a lot of effort to refine the previously schematic 
models with fluctuation inducing corrections and to em- 
ploy dynamical descriptions of the early non-equilibrium 
evolution pTHTij . 

Hydrodynamical simulations can take these fluctua- 
tions into account by generating an ensemble of runs 
each with a unique initial condition, so-called event by 
event simulations. This is in contrast to event averaged 
simulations where an ensemble of fluctuating initial con- 
ditions is generated, and then a single initial condition 
corresponding to this set's ensemble average is subject 
to evolution. Event by event modeling has proven to be 
essential for correctly describing all the details of the bulk 
behavior of heavy ion collisions [T5H25] . 

The two main models for the generation of hydro- 
dynamic initial conditions are the Glauber tl2j 26-28; 
and color glass condensate (CGC) models [25^151 . The 
Glauber model samples a Woods-Saxon nuclear density 
distribution for each nucleus. 

Color glass condensate models are ah initio calcula- 
tions motivated by the idea of gluon saturation of parton 



distribution functions at small momentum scales x. In 
CGC models the gluon distribution for each nucleon is 
computed and the nuclear collision is modeled as interac- 
tions between these coherent color flelds. Each of these 
models generates spatial fluctuations where the details 
depend on the assumptions used in the specific imple- 
mentation. Glauber fluctuations come from Monte-Carlo 
(MC) sampling the nuclear density distribution. CGC 
fluctuations arise similarly with additional contributions 
from the self interaction of the color fields. 

We present a method for generating a 2d decomposi- 
tion of fluctuations in the initial state energy density of a 
heavy-ion collision. We apply this framework to the en- 
sembles of initial states generated by UrQMD [55H37j . 
by an MC Glauber code of Qin [11], and by the MG- 
KLN code of Drescher and Nara |29l |30l |38| which is a 
based on CGC ideas. The events compared are generated 
for Au-f Au collisions at -y/s — 200 AGeV at two impact 
parameters & = 2, 7fm. We introduce summary statistics 
for the Fourier expansion and show how the radial infor- 
mation exposes clear differences between the fluctuations 
generated by UrQMD and the MC-KLN codes. 

Colliding nuclei are strongly Lorentz contracted along 
the beam axis, as such we neglect the longitudinal dimen- 
sion. A priori we expect fluctuations in both radial and 
azimuthal directions. Experiments make measurements 
in a 2d momentum space, and further analysis of the fl- 
nal state may reveal correlations with quantities derived 
from a fully 2d decomposition of the initial state. Recent 
work by |39j has shown that events can be constructed 
which have identical £2 , £3 but with dramatically different 
energy distributions leading to different flnal state flow 
coefficients. The structure of such events is not sensitive 
to an azimuthal only decomposition. 
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We seek a two dimensional Fourier expansion on a disk of radius tq. The orthogonality of Bessel functions of the 
first kind is the key to this decomposition 
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where Xa,n is the n*^ positive zero of Ja{x). For in- 
teger a, Ja{-X) = (-l)"Ja(x) and so [Ja+l{Xa,n)]'^ = 

[J|q.|+i(Aq,_„)]^, and Aa,„ = A_q,_„. It follows that the 
functions 
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form a complete orthonormal set on this disk (with the 
uniform measure r dr dO / tttq) . As such any well behaved 
function / on this disk admits the convergent expansion 
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in terms of these basis functions and a set of generalized 
Fourier coefficients „ € C, defined as 
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The basis functions (j)m,n are solutions to Bessel's equa- 
tion on the unit disk [40] . They are the eigenfunctions of 
the Laplacian with Dirichlet be and eigenvalue — (A^ „)^, 
such that 



The eigenvalues of the Laplacian A^„ „ define a character- 
istic inverse length scale for the associated Fourier com- 
ponent A„i n. Higher orders of m and n are associated 
with larger Xm,,n values corresponding to smaller spatial 
regions. A similar correspondence of length scales and 
Id Fourier coefficients has been pointed out in [41, 42J. 
A selection of these eigenvalues are plotted in Fig:[T] 
Other boundary conditions (e.g., Neumann) lead to sim- 
ilar countable sets of basis functions. 

Higher values of the angular indices ±to correspond to 
higher numbers of zero crossings in the angular compo- 
nents of basis functions, or "lumpiness" in rotation, while 
higher values of the radial indices n are associated with 
more roughness as one moves closer or farther from the 
center of mass. The first few basis functions are plotted 
in Fig:[2j The lumpy shape of the basis functions suggests 
that this representation will be very useful for character- 
izing the hot and cold spot structures in the initial state 
of a heavy ion reaction. 



A typical UrQMD event (after subtraction of the en- 
semble average) is shown alongside its decomposition and 
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FIG. 1: The eigenvalues A^ „ of the negative Dirichlet 
Laplacian on the disk are plotted as a function of |m| 
(angular) and n (radial). These can be interpreted as 
an inverse length scale for the coefficients Am,n- Note 
that in general \Xm+i,n - A,„.„| ^ \Xm,n+i - A,„^„|, 
although |Am,n| ^ 7'■(f^ + — 1/4) as n, m — )■ oo. 



the resulting Fourier coefficients Am,n in Fig:|3] Here we 
characterize the event in terms of coefficients m g [—8, 8] 
and n € [1,8], sufficient to capture a detailed image of the 
original initial state distribution. Higher order terms con- 
tribute very little additional information. The difference 
between the original distribution and the reconstructed 
distributions is dominated by numerical noise for coeffi- 
cients with order higher than m^^^ = 8, n^^x = 8. Let us 
now explore different ways to extract useful information 
from this decomposition. 

Using the orthogonality of the basis functions along 
with the Laplacian we can derive simple expressions for 
norms of the function / to be expanded in the frequency 
domain. The L2{f) norm, which measures the total mass 
of / is. 
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where (a, fe) = ^ a{r, 9)b{r, OydrdO is the inner 

product for functions on our disk. The Sobolev i?i(/) 
norm gives a measure of the total mass of the function 
as well as a measure of how 'wobbly' it is across the disk 
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FIG. 2: Plots of the first few real (top) and imaginary 
(bottom) components of 0m, n('',^')- The angular 
coefficient m G [0, 5] increases from left to right, the 
radial coefficient n € [1,3] increases from bottom to top. 



H,{f) :=((-^2^2+/)/,/)i/2 
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where ^ is a characteristic length scale introduced to 
maintain unit consistency (we use £ = Ifm). A varia- 
tion on the Sobolev norm gives the angular norm Mi{f) 
which quantifies angular gradients. 



1/2 



(7) 



We use these norms below to quantify properties of the 
events considered. 

To demonstrate the usefulness of our proposed method, 
we apply it to three example models: UrQMD, MC- 
Glauber and MC-KLN. Since this work does not aim at 
drawing definitive conclusions about initial state physics 
these models should be regarded as representative of the 
models currently in use. 

We consider a set of 100 events generated by each code. 
UrQMD includes Boltzmann hadronic transport before 
the hydro begins, while the MC-Glauber code includes 
simple streaming transport of the nucleons after inter- 
action; both of these will introduce added spatial fluc- 
tuation in the energy density. The Glauber code also 
includes KNO scaling of the multiplicity fluctuations per 



binary collision. For all models the initial condition was 
computed in a 200 point grid in the transverse plane at 
the center of the collision along the beam axis. To explore 
the centrality dependence of the analysis we consider 
events at two impact parameters b — 2,7fm. All events 
are generated for Au+Au collisions at -y/s^vjv = 200yl 
GeV. To study the fluctuations generated by these events 
the ensemble averaged event is computed for each model 
and subtracted from each event in the ensemble before 
applying the decomposition. 

In order to better understand the response of our pro- 
cess to fluctuations a series of progressively smoother en- 
sembles of UrQMD events were generated. The degree 
of fluctuation is controlled by populating these ensem- 
bles with events created by averaging N independent 
events together. See [l^ for more details on this pro- 
cess and its influence on the ellipticity and triangularity 
of the UrQMD initial state. As the number N of events 
over which we average increases, the result will increas- 
ingly resemble the ensemble average. We have examined 
events with TV = {1, 2, 5, 10, 25}. While only the iV = 1 
case represents the true UrQMD output, the diminishing 
amounts of fluctuation in the other successively smoother 
sets can be used to give perspective on (or "calibrate") 
our analysis. 

In Fig:|4]the ensemble averages of each norm {Hi, Mi, 
L2) are plotted for each set of events. Considering first 
the UrQMD results it is clear that each norm tends to 
zero as N increases as expected, confirming that each 
norm is quantifying fluctuations. Fluctuations are lower 
at the larger impact parameter b but the relative ordering 
of the models is preserved. For all norms we see good 
agreement between the N = 2 UrQMD events and the 
MC-Glauber code. 

The L2 norm computes the total amount of fluctua- 
tion, the Hi norm measures the spatial gradients in the 
fluctuations along with their total mass, the Mi norm 
measures angular gradients. The MC-KLN results show 
the largest Hi values while having Mi values compara- 
ble to UrQMD iV = 1, the MC-KLN L2 is somewhat 
less than UrQMD. The relative ordering of the Glauber 
and UrQMD events is the same across all norms. The 
UrQMD and Glauber events are rather similar in the 
scale and nature of their fluctuations, UrQMD produces 
slightly larger fluctuations. The MC-KLN model pro- 
duces fluctuations on a scale comparable to a hypotheti- 
cal UrQMD N = 3/2. However the large Hi and compa- 
rable Ml implies that these events exhibit larger radial 
gradients than the other models. This may be attributed 
to the very rapid spatial falloff of the gluon density near 
the edges of the nucleons. 

In Fig:[5] we plot the ensemble distributions of the 
norms. The width of the UrQMD distributions shrinks 
with increasing N and the modes of the distributions shift 
towards smaller values. The averaging procedure reduces 
fluctuations, and the norms are sensitive to these fluctua- 
tions. Across the range of norms and impact parameters 
the MC-Glauber results follow the UrQMD TV = 2 re- 
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FIG. 3: A typical UrQMD ensemble average subtracted event. Second: the reconstructed event generated by ^ 
after applying the decomposition. Bottom: The absolute values of the Fourier coefficients |Am,n| for this event. 
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FIG. 4: The ensemble averages of the Sobolev norm Hi (left) the angular norm Mi (center) and the L2 norm 
(right). The labels on the UrQMD glyphs give the number of averagings used to generate the ensemble members. 
The MC-Glauber and MC-KLN results are plotted with an artificial offset in b to permit easier comparison. 



suits very closely. The MC-KLN code produces events 
with a larger range of Hi values but with similar values 
to UrQMD N=l for the other angular norm and slightly 
reduced total mass (i2)- The MC-KLN events show a 
wider range of spatial and angular gradients than the 
Glauber based models. 

The ensemble average of the norms and their distribu- 
tions are instructive summaries of the fluctuations pro- 
duced by each model. In Fig:|6] we show the ensemble 
averages of the absolute values of the Fourier coefficients 
|^m,n| plotted against their eigenvalues Xm,n- These fig- 
ures should be thought of as generalized power-spectra; 
states with higher coefficients contribute more to the fluc- 
tuations than those with lower weights. The figures are 
plotted with coefficients of constant n (radial degree) 
joined by lines. For all models we see that the Fourier 
weights decrease rapidly with increasing angular degree 
|m|, along the lines, and more slowly with increasing ra- 
dial degree n at fixed m. These figures make it very clear 
that the different models produce quite different fluctu- 
ations. 

It is informative to contrast the largest few coefficients 
for each model. As shown by the norms the UrQMD 



and Glauber results are relatively similar — most of the 
weight is in coefficients with n = 1, 2 and |r7i| = 0, 4, cor- 
responding to events with rather more angular variation 
and less radial variation. The bulk of the mass of the 
distributions is at A ~ 5 — 10, implying that fiuctuations 
with large length scales dominate. 

The MC-KLN results stand in marked contrast. The 
largest weights are at n = 2,3,4 with A ~ 15, repre- 
senting basis states with more radial zero crossings and 
smaller characteristic sizes. We can conclude that the 
MC-KLN code generates fluctuations over smaller spa- 
tial scales with larger radial gradients compared to the 
Glauber based codes. These trends persist at both cen- 
tralities, although the absolute scales are diminished at 
larger b. This is reasonable given that there is a smaller 
active collision area at higher impact parameters. The 
2d Fourier decomposition is a straightforward method for 
capturing the differences in event structures in an intu- 
itive way. 

We have presented a new method for characterizing 
the fluctuations in the initial state of heavy ion collisions. 
The method is simple and general, it can be applied as 
easily to theoretical models as to the output of event gen- 
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FIG. 5: The distribution of the Hi (left), Mi (middle) and L2 (right) norms for UrQMD N = {1,2,5}, MC-KLN 
and MC-Glauber. The top row shows events at 5 = 2fm, the bottom row shows b — 7fm. 



erators. We have shown the ability of the norms we intro- 
duce to quantify the broad differences among the three 
models we considered. The radial information included 
by this process provides additional insights into the na- 
ture of fluctuations which are not readily attainable by 
considering quantities derived from angular decomposi- 
tions alone. 

We believe that these quantities and variants will prove 
to be useful tools for further understanding the differ- 
ences and similarities between initial condition models. 
They represent a reasonable compromise between too 
much information, as given by the full distributions of 
|>lm,n| above, and too little information as given by event 
shape measures which marginalize away the radial infor- 
mation. In future work we will examine how these struc- 
tures are passed through the hydrodynamical evolution 
to the hadronic final-state of the collision. Even if these 
quantities cannot be measured in detectors they provide 
a useful basis for apples-to-apples comparison of initial 
state models. 
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Appendix: Angular Decomposition 

To connect our analysis with commonly used quantities 
we quote in Table:[l]the values of 62 and £3 computed for 
events without ensemble average subtraction. 
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TABLE I: The ensemble average values of (£2) and (£3) 
for each of the models at impact parameters h ~ 2, 7fm. 
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FIG. 6: Comparison of |^m_„| (points) plotted against the characteristic length eigenvalue \ra,n at & = 2fm (top) 
and b = 7fni (bottom). The solid lines join points with constant n G [1, 8] the angular components included are in 
the range \m\ e [0, 8]. For a given n value the |m| distribution is monotonic, i.e., |m| = is always the leftmost point 

on the line and |m| = 8 is always the rightmost. 
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